#define LOGGING 1

#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <math.h>
#include "netcdf.h"
#include "HY1C_out.h"
#include "l12_proto.h"
#include "lonlat2pixline.h"
#include "l1_struc.h"
#include "l1_struc.h"
#include "l12_parms.h"
#include "h5io.h"
#include "config.h"
#include "hy1c_util.h"
#include "sqlite3_opt.h"

#define max(a,b)(a>b?a:b)
#define min(a,b)(a<b?a:b)

#define OUTPUT_REPORT 1

#define CHECK_AND_FREE_MEM(a) \
    do { \
        if (a != NULL) { \
            free(a); \
            a = NULL; \
        } \
    } while(0);

int32_t* year_hy1c_tp;
int32_t* day_hy1c_tp;

static int cal_day(int year, int days, int msec,
        int *month, int *day, int *hour,
        int *minutes)
{
    int i = 0;
    int days_month_orig[] = {0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
    int days_month_leap[] = {0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

    int leap = 0;
    if ((year % 100) == 0) {
        if ((year % 400) == 0) {
            leap = 1;
        }
    } else if ((year % 4) == 0) {
        leap = 1;
    } else {
        leap = 0;
    }

    int sum = 0;
    int *days_month = 0;
    if (leap) {
        days_month = days_month_leap;
    } else {
        days_month = days_month_orig;
    }

    for (i = 0; i < 12; i++) {
        sum += days_month[i];
        if ((days > sum) && (days <= sum + days_month[i + 1])) {
            *month = i + 1;
            *day = days - sum;
            //sprintf(date, "%d %s %d", year, month_str[i], days - sum);
            break;
        }
    }

    *hour = msec / 3600000;
    *minutes = (msec - (*hour)) / 1000;

    if ((*month <= 0) || (*month > 12)) {
        *month = 0;
    }

    if ((*day <= 0) || (*day > 31)) {
        *day = 0;
    }

    if ((*hour <= 0) || (*hour >= 24)) {
        *hour = 0;
    }

    if ((*minutes <= 0) || (*minutes >= 60)) {
        *minutes = 0;
    }

    return 0;
}

/**
 * @brief set logLevel and output file
 *
 * @param logfile logfile path, if NULL == logfile, output log to stdout
 *
 * @return 0 - success
 *         1 - failed
 */
int log_init(char *logfile)
{
    if (g_cfg.pLog)
        return 1;

    if (NULL == logfile) {
        g_cfg.log_level = LogLevel_Debug;
        g_cfg.pLog = stdout;
        return 0;
    }

    //////////////////////////////////////
    // ----------- LOGLEVEL -------------
    //////////////////////////////////////
    g_cfg.log_level = LogLevel_Debug;
    if ((g_cfg.pLog = fopen(logfile, "a")) == NULL) {
        fprintf(stderr, "Init log failed.\n");
        return 0;
    }

    return 1;
}

void log_uninit()
{
    fclose(g_cfg.pLog);
}

#define LOAD_L2_ITEM(a, b, c) \
    do { \
        if (b == NULL) \
        break; \
        a = b + c; \
    } while(0);

#define LOAD_L2_ITEM2(a, b, c, d) \
    do { \
        if (b == NULL) \
        break; \
        memcpy(a, b + c, d); \
    } while(0);

int loadl2_new(l2str *l2rec, l2rec_HY1C *data_hic, int iscan, int ipix, int type)
{
    init_l2_single(l2rec);

    int nwave = (type == 1) ? g_cfg.modis_nwave: g_cfg.hy1c_nwave;
    //l2rec->sensorID = g_cfg.sensor_id;
    l2rec->sensorID = HMODIST;
    l2rec->iscan = iscan;
    l2rec->ipix = ipix;
    l2rec->input = &(g_cfg.input);

    LOAD_L2_ITEM(l2rec->year,      data_hic->year,      iscan);
    LOAD_L2_ITEM(l2rec->day,       data_hic->day,       iscan);
    LOAD_L2_ITEM(l2rec->msec,      data_hic->msec,      iscan);
    LOAD_L2_ITEM(l2rec->invalid,   data_hic->invalid,   (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->Lr,        data_hic->Lr,        (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->TLg,       data_hic->TLg,       (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->tLf,       data_hic->tLf,       (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->Lw,        data_hic->Lw,        (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->La,        data_hic->La,        (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->Lt,        data_hic->Lt,        (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->Lt_theroy, data_hic->Lt_theroy, (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->taua,      data_hic->taua,      (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->tg_sol,    data_hic->tg_sol,    (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->tg_sen,    data_hic->tg_sen,    (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->t_sol,     data_hic->t_sol,     (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->t_sen,     data_hic->t_sen,     (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->gain,      data_hic->gain,      (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->aot_new,   data_hic->aot_new,   (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM(l2rec->Fo,        data_hic->Fo,        0);
    if (g_cfg.DPG_method == DPG_TRBMC) {
        LOAD_L2_ITEM(l2rec->field_AOT, data_hic->field_AOT, 0);
        LOAD_L2_ITEM(l2rec->field_Lw,  data_hic->field_Lw,  0);
        LOAD_L2_ITEM(l2rec->field_aer_modnum,  data_hic->field_aer_modnum,  0);
    }
    LOAD_L2_ITEM(l2rec->lat,       data_hic->lat,       (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->lon,       data_hic->lon,       (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->sena,      data_hic->sena,      (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->senz,      data_hic->senz,      (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->sola,      data_hic->sola,      (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->solz,      data_hic->solz,      (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->wv,        data_hic->wv,        (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->ws,        data_hic->ws,        (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->flags,     data_hic->flags,     (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->aermodmin, data_hic->aermodmin, (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->aermodmax, data_hic->aermodmax, (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->aerratio,  data_hic->aerratio,  (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->eps,       data_hic->eps,       (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->delphi,    data_hic->delphi,    (iscan * data_hic->npix + ipix));

    LOAD_L2_ITEM(l2rec->Tau_r,     data_hic->Tau_r,     0);
    LOAD_L2_ITEM(l2rec->rh,        data_hic->rh,        (iscan * data_hic->npix + ipix));
    LOAD_L2_ITEM(l2rec->pr,        data_hic->pr,        (iscan * data_hic->npix + ipix));

    static double esdist = 0;
    esdist = esdist_(l2rec->year,l2rec->day,l2rec->msec);
    l2rec->fsol = pow(1.0/esdist,2);

    if ((l2rec->sena == 0) || (l2rec->sola == 0)) {
        return 0;
    }

    *(l2rec->delphi) = *(l2rec->sena) - 180.0 - *(l2rec->sola);
    if (*(l2rec->delphi) < -180.)
        *(l2rec->delphi) += 360.0;
    else if (*(l2rec->delphi) > 180.0)
        *(l2rec->delphi) -= 360.0;

    return (0);
}

/* this method will copy the memory in fact */
int loadl2_new2(l2str *l2rec, l2rec_HY1C *data_hic, int iscan, int ipix, int type)
{
//    init_l2_single(l2rec);

    int nwave = (type == 1) ? g_cfg.modis_nwave: g_cfg.hy1c_nwave;
    //l2rec->sensorID = g_cfg.sensor_id;
    l2rec->sensorID = HMODIST;
    l2rec->iscan = iscan;
    l2rec->ipix = ipix;
    l2rec->input = &(g_cfg.input);

    LOAD_L2_ITEM2(l2rec->year,      data_hic->year,      iscan, sizeof(float));
    LOAD_L2_ITEM2(l2rec->day,       data_hic->day,       iscan, sizeof(float));
    LOAD_L2_ITEM2(l2rec->msec,      data_hic->msec,      iscan, sizeof(float));
    LOAD_L2_ITEM2(l2rec->lon,        data_hic->lon,       (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->lat,        data_hic->lat,       (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->solz,       data_hic->solz,      (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->sola,       data_hic->sola,      (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->senz,       data_hic->senz,      (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->sena,       data_hic->sena,      (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->delphi,     data_hic->delphi,    (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->Lt,        data_hic->Lt,        (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->aermodmin,  data_hic->aermodmin, (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->aermodmax,  data_hic->aermodmax, (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->aerratio,   data_hic->aerratio,  (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->eps,        data_hic->eps,       (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->taua,      data_hic->taua,      (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->TLg,       data_hic->TLg,       (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->La,        data_hic->La,        (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->Lw,        data_hic->Lw,        (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    //LOAD_L2_ITEM2(l2rec->invalid,   data_hic->invalid,   (iscan * nwave * data_hic->npix + ipix * nwave));
    LOAD_L2_ITEM2(l2rec->Lr,        data_hic->Lr,        (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->tLf,       data_hic->tLf,       (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->Lt_theroy, data_hic->Lt_theroy, (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->tg_sol,    data_hic->tg_sol,    (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->tg_sen,    data_hic->tg_sen,    (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->t_sol,     data_hic->t_sol,     (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->t_sen,     data_hic->t_sen,     (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->gain,      data_hic->gain,      (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    //LOAD_L2_ITEM2(l2rec->aot_new,   data_hic->aot_new,   (iscan * nwave * data_hic->npix + ipix * nwave), sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->Fo,        data_hic->Fo,        0, sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->wv,        data_hic->wv,        (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->ws,        data_hic->ws,        (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->Tau_r,     data_hic->Tau_r,     0, sizeof(float)*nwave);
    LOAD_L2_ITEM2(l2rec->rh,        data_hic->rh,        (iscan * data_hic->npix + ipix), sizeof(float));
    LOAD_L2_ITEM2(l2rec->pr,        data_hic->pr,        (iscan * data_hic->npix + ipix), sizeof(float));

    static double esdist = 0;
    esdist = esdist_(l2rec->year,l2rec->day,l2rec->msec);
    l2rec->fsol = pow(1.0/esdist,2);

    if ((l2rec->sena == 0) || (l2rec->sola == 0)) {
        return 0;
    }

    *(l2rec->delphi) = *(l2rec->sena) - 180.0 - *(l2rec->sola);
    if (*(l2rec->delphi) < -180.)
        *(l2rec->delphi) += 360.0;
    else if (*(l2rec->delphi) > 180.0)
        *(l2rec->delphi) -= 360.0;

    return (0);
}

int is_valid_point(l2str *l2, const int nwave)
{
    if (NULL == l2) {
        LOG_ERR("NULL is l2\n");
        return 0;
    }

    int i = 0;
    int is_invalid = 0;
    for (i = 0; i < nwave; i++) {
        if ((is_invalid |= l2->invalid[i]) != 0)
            printf("0x%X\n", l2->invalid[i]);
            break;
    }

    return !is_invalid;
}

void get_nearest_point(l2rec_HY1C *hy1c, l2rec_HY1C *modis, int32_t *hy1c_index, int32_t *modis_index)
{
    static float EARTH_RADIUS = 6378.137; // earch radius unit km.

    float rad(float d)
    {
        return d * M_PI / 180.0;
    }

    float get_distance(float lat1, float lng1, float lat2, float lng2)
    {
        float rad_lat1 = rad(lat1);
        float rad_lat2 = rad(lat2);
        float a = rad_lat1 - rad_lat2;
        float b = rad(lng1) - rad(lng2);

        float s = 2 * asin(sqrt(pow(sin(a / 2), 2) + cos(rad_lat1) * cos(rad_lat2) * pow(sin(b / 2), 2)));
        s = s * EARTH_RADIUS;
        return s;
    }

    float refer_lon = hy1c->field_AOT_lonlat[0];
    float refer_lat = hy1c->field_AOT_lonlat[1];
    int32_t nearest_index = 0;
    float nearest_distance = 0;
    float distance = 0;

    float lon = 0;
    float lat = 0;
    int i_h = 0;
    int i_m = 0;

    for (i_h = 0; i_h < hy1c->npix * hy1c->nscans; i_h++) {
        lon = hy1c->lon[i_h];
        lat = hy1c->lat[i_h];

        //distance = pow((lon - refer_lon), 2) + pow((lat - refer_lat), 2);
        distance = get_distance(refer_lat, refer_lon, lat, lon);
        if (i_h == 0) {
            nearest_distance = distance;
            nearest_index = i_h;
        } else {
            if (distance < nearest_distance) {
                nearest_distance = distance;
                nearest_index = i_h;
            }
        }
    }

    *hy1c_index = nearest_index;

    for (i_m = 0; i_m < modis->npix * modis->nscans; i_m++) {
        lon = modis->lon[i_m];
        lat = modis->lat[i_m];

        //distance = pow((lon - refer_lon), 2) + pow((lat - refer_lat), 2);
        distance = get_distance(refer_lat, refer_lon, lat, lon);
        if (i_m == 0) {
            nearest_distance = distance;
            nearest_index = i_m;
        } else {
            if (distance < nearest_distance) {
                nearest_distance = distance;
                nearest_index = i_m;
            }
        }
    }

    *modis_index = nearest_index;
}

void filter_gain(l2rec_HY1C *rec)
{
    float *Lt = rec->Lt;
    float *Lt_theroy = rec->Lt_theroy;
    float *gain = rec->gain;
    int32_t *invalid = rec->invalid;
    size_t nscan = rec->nscans;
    size_t npix = rec->npix;

    int32_t nwave = g_cfg.hy1c_nwave;
    size_t n = nscan * npix * nwave;

    size_t* pIndex[MAX_WAVE] = {0};
    float* pLt[MAX_WAVE] = {0};
    float* pLt_theroy[MAX_WAVE] = {0};

    for(int i = 0; i < nwave; i++) {
        pIndex[i] = (size_t*)malloc(nscan * npix * sizeof(size_t));
        pLt[i] = (float*)malloc(nscan * npix * sizeof(float));
        pLt_theroy[i] = (float*)malloc(nscan * npix * sizeof(float));
    }

    int32_t pIndex_pos[MAX_WAVE] = {0};
    int32_t Lt_valid_pos[MAX_WAVE] = {0};
    int32_t Lt_theroy_valid_pos[MAX_WAVE] = {0};

    int32_t i = 0, j = 0;
    for (i = 0; i < nwave; i++) {
        for (j = 0; i + j < (int32_t)n; j = j + nwave) {
            if (gain[i + j] == BAD_FLT)
                continue;

            //if (Lt[i + j] == BAD_FLT)
            //    continue;

            //if (Lt_theroy[i + j] == BAD_FLT)
            //    continue;

            pIndex[i][pIndex_pos[i]++] = i + j;
            pLt[i][Lt_valid_pos[i]++] = Lt[i + j];
            pLt_theroy[i][Lt_theroy_valid_pos[i]++] = Lt_theroy[i + j];
        }
    }

    for (int i = 0; i < nwave; i++) {
        float stdev_Lt = stdev_float(pLt[i], Lt_valid_pos[i]);
        float mid_Lt = getmid_float(pLt[i], Lt_valid_pos[i]);

        float stdev_Lt_theroy = stdev_float(pLt_theroy[i], Lt_theroy_valid_pos[i]);
        float mid_Lt_theroy = getmid_float(pLt_theroy[i], Lt_theroy_valid_pos[i]);

        for (int j = 0; j < Lt_valid_pos[i]; j++) {
            if (fabs(pLt[i][j] - mid_Lt) > stdev_Lt / 2) {
                gain[pIndex[i][j]] = BAD_FLT;
                invalid[pIndex[i][j]] |= 0x40000000;
                LOG_DBG("[%d][%04d][%04d] invalid Lt %f mid %f stdev %f\n", i, j / npix, j % npix,  pLt[i][j], mid_Lt, stdev_Lt);
                continue;
            }

            if (fabs(pLt_theroy[i][j] - mid_Lt_theroy) > stdev_Lt_theroy) {
                gain[pIndex[i][j]] = BAD_FLT;
                invalid[pIndex[i][j]] |= 0x40000000;
                LOG_DBG("[%d][%04d][%04d] invalid Lt_theroy %f mid %f stdev %f\n", i, j / npix, j % npix,  pLt_theroy[i][j], mid_Lt_theroy, stdev_Lt_theroy);
                continue;
            }
        }
    }

    for(int i = 0; i < nwave; i++) {
        free(pIndex[i]);
        free(pLt[i]);
        free(pLt_theroy[i]);
    }

    return;
}

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <dirent.h>
#include <unistd.h>

extern l2rec_HY1C* read_H1C_l2_hdf5(char* filename, int type);

/**
 * @brief the entry for calculate the gain
 *
 * @param result[] result of gain
 * @param length gain length
 *
 * @return 0 - success
 *         1 - fail
 */
int coeff_calc(float result[], size_t length)
{
    int ret = 0;

    static int first = 1;
    static l2rec_HY1C *l2recHY1C_rid = NULL;
    static l2rec_HY1C *l2recMODIS_rid = NULL;

    const char * inHY1C_F = g_cfg.hy1c_f;
    const char * inMODIS_F = g_cfg.modis_f;
    const char * configFile = g_cfg.config_f;
    char * outGain = g_cfg.out_gain;
    char** aermodules = g_cfg.aermodules;
    int naermodules = g_cfg.naermodules;

    float *wave = g_cfg.wave_hy1c;
    int32_t opoints = 0; // all the points enjoy the cal.
    int32_t nwave = g_cfg.hy1c_nwave;

    if (first) {
        LOG_INFO("HY1C input file is %s\n",  inHY1C_F);
        LOG_INFO("MODIS input file is %s\n", inMODIS_F);
        LOG_INFO("configFile path is %s\n",  configFile);
        LOG_INFO("output file path is %s\n", outGain);

        LOG_DBG("Begin to read HY1C data.\n");
        l2recHY1C_rid = read_H1C_l2_hdf5(inHY1C_F, 0);

        LOG_DBG("Begin to read MODIS data.\n");
        l2recMODIS_rid = read_H1C_l2_hdf5(inMODIS_F, 1);
    }

    l2str l2rec[2] = {0};
    l2str* l2recHY1C = &l2rec[0];
    l2str* l2recMODIS = &l2rec[1];

    int32_t iscan_M = 0;
    int32_t iscan_H = 0;
    int32_t ipix_M = 0;
    int32_t ipix_H = 0;

    int32_t sscan_M = 0;
    int32_t spix_M = 0;
    int32_t nscan_M = l2recMODIS_rid->nscans;
    int32_t nscan_H = l2recHY1C_rid->nscans;
    int32_t npix_M = l2recMODIS_rid->npix;
    int32_t npix_H = l2recHY1C_rid->npix;
    int32_t hy1c_coordinate = 0;
    int32_t valid_point_count  = 0;
    int32_t ngain = nscan_H * npix_H * nwave;

    /* these below variables only used for TRBMC */
    int32_t trbmc_iscan_M = 0;
    int32_t trbmc_ipix_M = 0;

    if (g_cfg.DPG_method == DPG_TRBMC) {
        int32_t hy1c_index = 0;
        int32_t modis_index = 0;
        get_nearest_point(l2recHY1C_rid, l2recMODIS_rid, &hy1c_index, &modis_index);

        /*
         * calc all the points in the rectangular region centered at this
         * particular point within 2n. if n = 5, the points are 11*11.
         */

        /*
         * for TRBMC no need to get the hy1c coordinate according to modis', so
         * iterate hy1c directory.
         */

        int n = g_cfg.trbmc_region;

        iscan_H = hy1c_index / npix_H;
        ipix_H = hy1c_index % npix_H;

        trbmc_iscan_M = modis_index / npix_M;
        trbmc_ipix_M = modis_index % npix_M;

        LOG_INFO("TRBMC matched : ISCAN_H %u IPIX_H %u ISCAN_M %u IPIX_M %u\n",
                iscan_H, ipix_H, trbmc_iscan_M, trbmc_ipix_M);

        int32_t scan_start_hy1c = (iscan_H - n) > 0 ? (iscan_H - n) : 0;
        int32_t pixel_start_hy1c = (ipix_H - n) > 0 ? (ipix_H - n) : 0;

        int32_t scan_end_hy1c = (iscan_H + n + 1) < nscan_H ? (iscan_H + n + 1) : nscan_H;
        int32_t pixel_end_hy1c = (ipix_H + n + 1) < npix_H ? (ipix_H + n + 1) : npix_H;

        sscan_M = scan_start_hy1c;
        spix_M = pixel_start_hy1c;

        nscan_M = scan_end_hy1c;
        npix_M = pixel_end_hy1c;
    }

#ifdef OUTPUT_REPORT
    init_db(NULL);
#endif

    /* store iscan_M ipix_M for trbmc and restore them when finish the loop */
    int32_t store_iscan_M = 0;
    int32_t store_ipix_M = 0;

    l2str l2rec2;
    alloc_l2_single(8, &l2rec2);

    for (iscan_M = sscan_M; iscan_M < nscan_M; iscan_M++) {
        for (ipix_M = spix_M; ipix_M < npix_M; ipix_M++) {
            if (g_cfg.DPG_method != DPG_TRBMC) {
                hy1c_coordinate = l2recHY1C_rid->modis_matched_flag[iscan_M * npix_M + ipix_M];
                if (hy1c_coordinate == -1) {
                    /* There is no mapped HY1C data for this point */
                    continue;
                }

                if ((iscan_H = hy1c_coordinate / 100000) > nscan_H) {
                    LOG_ERR("scan %d is out range %d\n", iscan_H, nscan_H);
                    continue;
                }

                if ((ipix_H = hy1c_coordinate % 100000) > npix_H) {
                    LOG_ERR("pix %d is out range %d\n", ipix_H, npix_H);
                    continue;
                }

                if (iscan_H > 20000) {
                    LOG_ERR("scan %d is out range %d\n", iscan_H, nscan_H);
                    continue;
                }

                if (ipix_H > 99999) {
                    LOG_ERR("pix %d is out range %d\n", ipix_H, npix_H);
                    continue;
                }
            } else {
                store_iscan_M = iscan_M;
                store_ipix_M = ipix_M;

                iscan_H = iscan_M;
                ipix_H = ipix_M;

                iscan_M = trbmc_iscan_M;
                ipix_M = trbmc_ipix_M;
                LOG_INFO("iscan_M [%d] ipix_M [%d] iscan_H [%d] ipix_H [%d]\n", iscan_M, ipix_M, iscan_H, ipix_H);
            }

            if (loadl2_new(l2recHY1C, l2recHY1C_rid, iscan_H, ipix_H, 0) != 0) {
                LOG_ERR("Loadl2 failed\n");
                /* restore TRBMC */
                if (g_cfg.DPG_method == DPG_TRBMC) {
                    iscan_M = store_iscan_M;
                    ipix_M = store_ipix_M;
                }
                continue;
            }

            if (loadl2_new(l2recMODIS, l2recMODIS_rid, iscan_M, ipix_M, 1) != 0) {
                LOG_ERR("Loadl2 failed\n");
                /* restore TRBMC */
                if (g_cfg.DPG_method == DPG_TRBMC) {
                    iscan_M = store_iscan_M;
                    ipix_M = store_ipix_M;
                }
                continue;
            }

            /* restore TRBMC */
            if (g_cfg.DPG_method == DPG_TRBMC) {
                iscan_M = store_iscan_M;
                ipix_M = store_ipix_M;
            }

            opoints++;

#ifndef FD_DEBUG
            /* only concern about the valid gain points for the second calc. */
            if (!first && l2recHY1C->gain[nwave - 1] == BAD_FLT) {
                LOG_DBG("HY1C[%04d][%04d] MODIS[%04d][%04d] HY1C gain is not valid\n",
                        iscan_H, ipix_H, iscan_M, ipix_M);
                continue;
            }

            if (!is_valid_point(l2recHY1C, g_cfg.hy1c_nwave)) {
                LOG_DBG("HY1C[%04d][%04d] MODIS[%04d][%04d] HY1C is not a valid point\n",
                        iscan_H, ipix_H, iscan_M, ipix_M);
                continue;
            }

            if (!is_valid_point(l2recMODIS, g_cfg.modis_nwave)) {
                LOG_DBG("HY1C[%04d][%04d] MODIS[%04d][%04d] MODIS is not a valid point\n",
                        iscan_H, ipix_H, iscan_M, ipix_M);
                continue;
            }

            if ((*(l2recMODIS->flags) & 0x21b) > 0) {
                LOG_DBG("HY1C[%04d][%04d] MODIS[%04d][%04d] MODIS flag is not valid\n",
                        iscan_H, ipix_H, iscan_M, ipix_M);
                continue;
            }
#endif
            if (first) {
                ret = HY1C_cal_La(l2recHY1C, l2recMODIS, wave, nwave,
                                  aermodules, naermodules, configFile);
            } else {
                loadl2_new2(&l2rec2, l2recHY1C_rid, iscan_H, ipix_H, 0);
                ret = HY1C_cal_La2(result, length, &l2rec2, l2recMODIS, NULL);
                memcpy(l2recHY1C->La, l2rec2.La, sizeof(float) * nwave);
            }

            int invalid_la = 0;
            for (int iw = 0; iw < nwave; iw++)
            {
                if (((l2recHY1C->La)[iw] < 0) || ((l2recHY1C->La)[iw] > 10)) {
                    l2recHY1C->gain[iw] = BAD_FLT;
                    invalid_la = 1;
                    break;;
                }
            }

            if (invalid_la)
                continue;

            valid_point_count++;

            for (int i = 0; i < nwave; i++) {
                l2recHY1C->gain[i] = BAD_FLT;
            }

            ret = HY1C_cal_gain2(l2recHY1C, l2recMODIS, wave, nwave);

#ifdef OUTPUT_REPORT
            insert_report(l2recHY1C, l2recMODIS, wave, nwave);
#endif
        }
    }

    filter_gain(l2recHY1C_rid);

#ifdef OUTPUT_REPORT
    LOG_INFO("Begin to generate gain\n");
    insert_gain_report(l2recHY1C_rid, l2recMODIS_rid, wave, nwave);
    generate_report();
    close_db();
#endif

    // get the average per wave
    // for CZI, there are 4 waves totally and 2 cameras, should split the gain
    // data from left to right. gain[0-3] store the gain for the camera on the
    // left and gain[4-7] store the gain for the camera on the right.
    // index 0  1   2   3   4   5   6   7
    // gain  l  l   l   l   r   r   r   r
    float *gain = l2recHY1C_rid->gain;
    float gain_sum[MAX_WAVE] = {0};
    int32_t gain_count[MAX_WAVE] = {0};

    if (g_cfg.sensor_id == SENSOR_CZI) {
        int32_t *side_flag = l2recHY1C_rid->side_flag;
        int32_t start_l = side_flag[0];
        int32_t end_l = side_flag[1];
        int32_t start_r = side_flag[2];
        int32_t end_r = side_flag[3];

        int32_t i = 0, j = 0;
        for(i = 0; i < nwave; i++) {
            for(j = 0; i + j < ngain; j = j + nwave) {
                if (gain[i + j] == BAD_FLT)
                    continue;

                int point_index = j / nwave % npix_H;
                if ((point_index >= start_l) && (point_index <= end_l)) {
                    gain_sum[i] += gain[i + j];
                    gain_count[i]++;
                }

                if ((point_index >= start_r) && (point_index <= end_r)) {
                    gain_sum[i + 4] += gain[i + j];
                    gain_count[i + 4]++;
                }

            }
        }
    } else {
        int32_t i = 0, j = 0;
        for(i = 0; i < nwave; i++) {
            for(j = 0; i + j < ngain; j = j + nwave) {
                if (gain[i + j] == BAD_FLT)
                    continue;

                gain_sum[i] += gain[i + j];
                gain_count[i]++;
            }
        }
    }

#ifndef FD_DEBUG
    if (valid_point_count <= 0) {
        LOG_ERR("no valid point for cal.\n");
        return 1;
    }
#endif

    int year = *(l2recHY1C_rid->year);
    int month = 0;
    int day = 0;
    int hour = 0;
    int minutes = 0;

    cal_day(year, *(l2recHY1C_rid->day), *(l2recHY1C_rid->msec), &month, &day, &hour, &minutes);
    //cal_day(year, *(l2recHY1C_rid->day), 0, &month, &day, &hour, &minutes);

    FILE * fp = NULL;

    char output[256] = {0};
    if (first && g_cfg.again.keep_original_gain) {
        sprintf(output, "%s0", outGain);
    } else {
        sprintf(output, "%s", outGain);
    }

    if ((fp = fopen(output, "w")) == NULL) {
        printf("%s\n","can not open file");
        return 1;
    }

    fprintf(fp, "%s\n", "/begin_header");
    fprintf(fp, "%s\n", "/Satellite=HY-1C");
    fprintf(fp, "%s%s\n", "/Sensor=", g_cfg.sensor);
    fprintf(fp, "/Date=%04d%02d%02d\n", year, month, day);
    fprintf(fp, "/Time=%02d%02d\n", hour, minutes);
    fprintf(fp, "%s%d\n", "/All Points=", opoints);
    fprintf(fp, "%s%d\n", "/Valid Points=", valid_point_count);
    fprintf(fp, "%s%0.2f\%\n", "/Valid Points/All Points=",
            (float)valid_point_count/opoints * 100);
    fprintf(fp, "%s%s\n", "/Calibration_Method=", g_cfg.DPG_method_str);
    fprintf(fp, "%s%d\n", "/BandNumber=", g_cfg.hy1c_nwave);
    fprintf(fp, "%s\n", "/fields=Band,Coefficient");
    fprintf(fp, "%s\n", "/line_number=8");
    fprintf(fp, "%s\n", "/column_number=2");
    fprintf(fp, "%s\n", "/end_header");

    for(int i = 0; i < length; i++) {
        result[i] = gain_sum[i] / gain_count[i];
        if (!isnormal(result[i]))
            result[i] = BAD_FLT;

        fprintf(fp, "%d\t%0.1f\n", i+1, result[i]);
    }

    fclose(fp);

    first = 0;
    return 1;
}

char* get_version(char *arg)
{
    static char version[128] = {0};
    sprintf(version, "\n\n%s\nVersion: 1.3.3\nBuild: %s %s\n\n", arg, __DATE__, __TIME__);
    return version;
}

#include <libgen.h>
int main(int argc, char* argv[])
{
    log_init(NULL);
    LOG_INFO(get_version(basename(argv[0])));

    if (init_cfg(argc, argv) == 1) {
        log_uninit();
        return 0;
    }

    float result[8]; // TODO hardcode to consider with CZI
    memset(result, 0x00, sizeof(result));

    int ret = coeff_calc(result, sizeof(result)/sizeof(result[0]));

    if (g_cfg.again.calc_again) {
        ret = coeff_calc(result, sizeof(result)/sizeof(result[0]));
    }

    LOG_INFO("Finished -------------------------\n");
    log_uninit();
    return ret;
}
